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Abstract 

We performed first-principles molecular dynamics calculations for lithium using the projector 

O 

' augmented waves method and the generalized gradient approximation as exchange-correlation en- 

ergy. The melting curve of lithium was computed using the Z-method technique for pressures up to 
30 GPa, which agrees well with the experimental and two-phase simulated results. The change of 
^ ' the melting line slope from positive to negative was predicted by the characteristic shape inversion 

■ of the Z curve at about 8.2 GPa. Through analyzing the static properties, we conclude that no 

O . liquid-liquid phase transition accompanies the occurrence of the melting line maximum, which is 

caused by the higher compressibility of the liquid phase compared to the solid phase. In addition, 
^ ■ we systematically studied the dynamic and optical properties of lithium near melting curve at crit- 

^ : 

0^ ■ ical superheating and melting temperatures. It was suggested that spectra difference at critical 

(N 

I superheating and melting temperature may be able to diagnose the homogeneous melting. 

d ' 

O '. PACS numbers: 71.27.+a, 71.15.Mb, 71.20.-b, 63.20.dk 
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I. INTRODUCTION 



The simple alkali metals like Li and Na have drawn extensive attentions recently due 
to their enigmatic melting behavior at high pressure. Intensive theoretical and experimen- 
tal studies have demonstrated that they undergo a sequence of symmetry-breaking struc- 
tural phase transitions and present rather complex crystalline phases under compression3i~-. 
These structural changes are accompanied by a variety of intriguing phenomena, among 
which the most striking is the anomalous melting feature at high pressure. For sodium, 
its melting curve has been measured up to 130 GPa^ and subsequently reproduced by 
first-principles calculations based on either molecular dynamics or the usual Lindemann 
criterion^"-. These studies have revealed the unusual melting behavior of Na, i.e., the ex- 
istence of multiple maxima. For lithium, on the other hand, because of the difficulties in 
containing the sample under high pressure, the knowledge of its melting curve has long 
been confined to be less than 8 GPa^^ until a recent differential thermal analysis (DTA) 
measurement's^ which extended the melting line of Li up to 15 GPa and reported a maxi- 
mum at about 10 GPa. Theoretically, a Lindermann model curve of Li was calculated^ to 
give an obvious discontinuity near the bcc-fcc-liquid triple point, which was not supported 
by experimental data^. More directly, Tamblyn et al.^^ and Hernadez et alM have from 
first-principles molecular dynamics (FPMD) simulations determined the melting tempera- 
ture of Li over a broad pressure range. They both observed the negative slope of the melting 
curve. In addition, a new phase in liquid lithium with sp'^ bonded tetrahedral local order at 
pressures above 150 GPa was also predictecU^. A possible link between the maximum in the 
melting line and liquid-liquid phase transition (LLPT) was suggested for some systems, such 
as pi3,i4 g^j^^j Cs^^. As for Na, first-principles simulations have revealed that the maximum 
in the melting line occurs without any accompanying LLPT, and higher compressibility of 
the liquid phase than the solid phase causes the change of melting line slope from positive to 
negative^. The cause of the anomalous melting behavior of Li still remains unclear. Though 
previous FPMD studies have found bcc-like to fcc-like structural transition in liquid Li^^, 
the temperatures are well above the melting temperature and thus it could not be concluded 
that the similar structural transition occurs along the melting line. It is desirable to explore 
whether LLPT exists along the melting line and reveal what contributes to the melting curve 
maximum of Li. 
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There are two common strategies for FPMD calculation of the melting curve, i.e., the two- 
phase method^ii^ corresponding to a heterogeneous melting and the one-phase approaches 
which involves a homogeneous melting. The first one can give accurate results at the cost of 
computational demanding, while the second usually causes superheating effect, though needs 
much less atoms. A useful alternative has been proposed by Belonoshko et al.—, which is 
referred to as Z-method. It combines the advantages of both techniques. When the system 
reaches the critical superheating Tis, it could not be heated before transforming into a liquid 
structure. If let it evolve naturally as in Z-method, the temperature just drops down to 
the melting temperature T^. Along this line, recently, growing interests are concentrated 
on the mechanism of homogeneous melting through characterizing the crystal properties 
at Tis—'—. However, no consensus could be reached. The dynamic and optical properties 
of crystal at Tis, which are deemed useful for insight into the mechanism of melting, are 
still scarcely presented in the literature. Besides, the optical properties should be different 
for the superheating solid and the disordered liquid phase, and thus they could be able to 
diagnose the homogeneous melting. The similar idea has been suggested for diagnosing the 
shock melting of Al^^. 

Inspired by the above-mentioned facts, in this paper we calculate the melting curve of 
Li up to 30 GPa using the Z-method implemented by FPMD simulations. We show that 
the shape of Z curve inverse with pressure increased to ~8.2 GPa, which predicts the pres- 
ence of maximum in the melting curve, in good agreement with experimental measurement. 
Through examining the static properties of Li at Tis and Tm, we conclude that no LLPT 
occurs and liquid phase is more compressible than solid phase, which may be the cause to 
the melting line maximum. Besides, the dynamical and optical properties at and T^g 
are studied. It is found that solid and liquid spectra show marked difference. This could be 
an efficient way of diagnosing the phase transition during the homogeneous melting. In the 
next section, the methods used in homogeneous melting simulation and optical properties 
calculations are described. In Sec. Ill, the calculated results are discussed and compared 
with experimental data. Finally, we close our paper with a summary of our main results. 
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II. METHOD 



We have performed FPMD simulations to determine the melting curve using Z-method, 
which has been proven successful to predict melting temperature in several systems, such 
as A&; H2I, MgO^^ and so on. The idea is to perform FPMD in the microscopic ensemble 
(NVE) on a single solid system at different initial K (kinetic energy). A realistic Tis can 
be reached without any external intervention on the dynamics of the melting process. On 
further increasing K slightly, Tis will drop naturally to the melting temperature Tm at 
the pressure fixed by the chosen density. By performing long microscopic simulations at 
different cell volumes, it is possible to obtain points (P, T^) directly on the melting curve. 
The method is as straightforward as the two-phase method, and it requires half as many 
atoms in the simulation cell. However, it still requires a large amount of simulation steps to 
achieve complete melting curve. 

We performed Z-method simulations of Li melting with ab initio Simulation Package 
(VASP)^ for bcc structure (for eight densities). The all-electron projector augmented wave 
(PAW) method^i^ was adopted, retaining only the 2s electron in the valence, and the 
exchange-correlation energy was described employing the generalized gradient approxima- 
tion (GGA) of Perdew-Burke-Ernzerhof (PEE) formalism^^, as implemented in VASP. We 
used plane-wave cutoff of 150 eV, and the Brillouin zone was sampled only with the F point. 
For each density the system was simulated for 10000-20000 steps with the time step of 2.0 
fs, where the time scale lies between 20 to 40 ps for different initial K in order to construct 
an isochoric curve P versus T. In this study, although the applied pressure range is up to 
30 GPa, only the bcc structure was used in the MD calculations. A 256-atom bcc supercell 
was constructed. 

Following the FPMD simulations, a total of ten configurations are selected from an equi- 
librated (in an average sense) portion of the molecular dynamics run, typically sampling the 
final picosecond of evolution. For each of these configurations, the electrical conductivity is 
calculated using the Kubo-Greenwood formulai^^. The Kubo-Greenwood formulation gives 
the real part of the electrical conductivity as a function of frequency u, 
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Table I: Melting temperature and critical superheating temperature Tig for Li at different 
pressures. 

Melting Superheating 

Pressure 

temperature temperature 



(GPa) 


(K) 


(K) 


1.31±0.03 


495.36±9.27 


578.16±13.15 


3.18±0.03 


510.29±9.73 


613.38±12.46 


5.35±0.02 


521.76±10.13 


608.36±12.11 


8.21±0.02 


541.15±10.68 


627.1±13.16 


10.14±0.02 


538.78±10.34 


626.82±12.86 


15.69±0.02 


502.62±9.21 


591.42±12.02 


19.60±0.04 


489.32±7.92 


578.78±10.99 


28.58±0.04 


411.56±7.92 


480.20±8.93 



where /(ei,k) describes the occupation of the ith band, with the corresponding energy ej^k 
and the wave function ^j^k at k, and w (k) is the k-point weighting factor. 

III. RESULTS AND DISCUSSION 
A. Melting curve 

By performing the microscopic Z-method simulations, the system evolve freely without 
any temperature control, and in each case, after reaching TJ^, the isochore line drops to a 
point (P, Tjn) that should fall precisely along the melting curve. The isochore plots for each 
density are shown in Fig. 1, from which the melting points (the ones marked with square 
on each plot) and critical superheating points are extracted and shown in Table I with their 
respective error estimations. 

As can be seen in the plots in Fig. 1, these eight isochores naturally fall into two categories 
according to their characteristic shapes. The first three form a "Z" shape, while the last five 
form an inverse "Z" shape. The fourth plot should be noticed that the upper and lower 
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Figure 1: Isochoric curves for p = 0.58, 0.65, 0.71, 0.78, 0.82, 0.93, 1.0, and 1.13 g/cm^. The 
respective melting points are marked as squares in the different plots. Solid lines are just a guide 
for eye. 



cap of "Z" are nearly overlapped. These featured isochore characteristics just predict the 
anomalous melting behavior of lithium. When the volume is fixed at values corresponding 
to the isochores with "Z" shape, the pressure of the liquid phase is higher than that of the 
solid phase at the melting temperature. This implies that in condition of constant pressure, 
the liquid phase would have a larger volume, and thus the pressure derivative of the melting 
line in this region should be positive according to the Clausius-Clapeyron equation 



dP ~ ""AH' ^ ^ 

where is the melting temperature, P is the pressure, A V = Vi — Vg is the difference of 
molar volumes, and A H — Hi — Hg is the difference of molar enthalpies. On the other hand, 
in the region of the isochores with inverse "Z", the case is totally opposite. That is, the sohd 
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Figure 2: Pressure difference between tlie liquid [Pi) and solid (P^) phases as a function of volume. 
The corresponding pressures are indicated. 

phase would have a larger volume, consistent with the negative slope of the melting curve. 
In the case that the "Z" shape almost merges into one line as shown in the fourth plot in 
Fig. 1, the melting temperature would reach its maximum. We further show the pressure 
difference between the liquid (Pi) and solid (P^) phases as a function of volume at melting 
point in Fig. 2, which are obtained from the isochore simulations. It is indicated that the 
melting curve maximum locates at about 8.2 GPa, where Pi equals to P^. This agrees well 
with the recent experimental measurement^ of a negative slope above ~lo GPa. 

The melting curve along with the estimated errors is shown in Fig. 3, in comparison with 
experimental data and data from two-phase FPMD simulations. The maximum melting 
temperature of ~540 K is determined to be in the pressure range from 8 to 10 GPa. There is 
reasonable agreement between the Z-method melting and experimental results from Boehler 
et al? and Luedemann et alA, where the discrepancy is below 20 K. Especially, they reports 
a melting temperature of 508 K at 3 GPa, while we obtain 510 K at 3.18 GPa. And also 
unlike other single-phase simulations of Li, the present Z-method produces melting line 
very close to the two-phase simulation results up to 30 GPa we considered. The two-phase 
simulations have shown an interesting feature that bcc is more stable than fee structure 
close to meltingi^. Therefore, it should be reasonable to study the properties of lithium at 
melting and superheating temperatures just based on the simulation results of bcc structure 
in the following. 

First-principles simulations have confirmed that the anomalous melting behavior of Na 
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Figure 3: (Color online) Melting curve for lithium up to 30 GPa. The squares represent the 
individual points obtained with the Z-method. The diamonds are some experimental points from 
Luedemann et al. (Ref. 11) up to 8 GPa, crosses are experimental points from Bochler et al. (Ref. 
12) up to 3 GPa, and circles are recent experimental points from Lazicki et al. (Ref. 13) up to 15 
GPa. The up-triangles and down-triangles are the theoretical points of bcc and fee with two-phase 
simulations, respectively (Ref. 15). 

is attributed to the high compressibility of liquid phase than soUd phase, and the maximum 
of the melting line occurs without any accompanying first-order LLPT. Through analyzing 
the pair correlation function (PCF) 5'(r), we find that these conclusions also hold for Li. 
As we all know, PCF is usually used to examine atomic configurations, defined as pg (r) = 

( XI ^ l'" + ~ ^i) /' with r the interatomic distance, the number of atoms, p the 
density N/V , and and Vj the positions of atoms i and j, respectively. In Fig. 4(a) density 
dependence of PCF at and Tis are presented with r scaled by tq, where rQ={N /VY^'^ . 
As seen in the figure, the shape of the PCF does not change up to 1.13 g/cm^ (~30 GPa), 
which indicates that the compression in this pressure range is uniform with the local structue 
unchanged. However, the structure differences between Tis and Tm are noticeable especially 
for large values of r. For example, the peak structure at around r/r^ = 2.8 at Ti^ is large 
blurred at T^. Furthermore, the coordination at and Ti^ for different densities are given 
from integrating g{r) within a sphere of radius i?, as shown in Fig. 4(b). In particular, 
by integrating g (r) to the first minimum, the average coordination numbers Cnn can be 
obtained, as indicated in the Fig. 4(b). The first minimum of g (r) at is at r/ro~1.5, 
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Figure 4: (Color online) Density dependence of (a) pair correlation function g (r) and (b) coordi- 
nation number at (solid line) and Ti^ (short dashed line). The vertical dotted lines indicate the 
first minimum of g (r). To facilitate the comparison between different densities, r is scaled by tq. 

which is indicated by the dotted hne. It is obvious that Cat a? is nearly constant at 14 for 
these densities, similar to that of a bcc crystal. Thus our simulations give that no LLPT 
occurs along the melting line of Li up 30 GPa. It is consistent with the experimental results^ 
that no signature of bcc-fcc transition along the melting line was observed up to 15 GPa, 
though liquid structural transformation from bcc to fee in pure liquid phase above melting 
temperature is a natural case^^. In addition, as can be seen from Fig. 4(b), beyond the first 
coordination shell, the coordination number for the liquid phase increases compared to the 
solid phase, which implies that the liquid phase is more compressible than the solid phase, 
and thus may lead to the maximum of the melting line. 

Besides the above-presented clarification of the structure, we further examine the mean 
square displacement (MSD) at Tm and Tis in order to confirm that we indeed have a solid 
behavior at T/j, and a liquid behavior at for each of the volumes considered. Figure 5 
shows the MSD up to 3 ps for P ~ 10 GPa at Tis = 627 K and = 539 K. For purpose 
of comparison, we also include the MSD at T = 540 K (solid) and T = 571 K (liquid). It 
is obvious to see that the MSD at Tm increases linearly at long time, in a similar way as 
in pure liquid at T = 571K, while at Tis the displacement reaches a constant nearly, which 
suggests that atoms are in a solid structure and could not diffuse away from their equilibrium 
positions (compare with MSD at T = 540 K). 
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Figure 5: (Color online) MSD of Li atoms for P ~ 10 GPa at T = 540 K (dashed line), T;^ = 627 
K (dotted line), Tm = 539 K (thin solid line), and T = 571 K (thick solid line). Inset shows the 
same curves up to 0.3 ps. 

The diffusion coefficient D, estimated using the MSD up to to = 0.5 ps by 



is shown in Fig. 6 for every pressure point in our calculated melting curve. Here it is verified 
that there is a remarkable difference in the atomic diffusion between at Tm and at Tis in 
the whole pressure range considered. As expected, the diffusion decreases with pressure 
generally. Especially, the diffusion coefficient at 1.31 GPa and = 495 K is about 0.69 
A^/ps, which is in excellent agreement with the experimental value of 0.69 ± 0.09 A^/ps at 
GPa and 470 K^. 



B. Dynamic conductivity and optical properties 

The linear optical conductivities of superheated and melted Li at different densities (0.58, 
0.78 and 1.13 g/cm^) calculated using Eq. (1) are shown in Fig. 7, which are averaged over 
ten snapshots selected during the course of FPMD simulations. For identical density, the 
spectra at Tis and Tm show marked differences. At Tis, there are some structural peaks, 
while the dips fill in and only leave a shoulder at Tm- Of note is the fact that our predicted 
difference between superheated solid and melting liquid is even more pronounced at lower 
densities. For example, at density of 0.58 and 0.78 g/cm^, there are two prominent peaks, 
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Figure 6: (Color online) Diffusion coefficient for Li atoms compared for each pressure at melting 
point (circles) and superheating point (squares). 




Photon energy (eV) 

Figure 7: (Color online) The lithium optical conductivity versus to at (a) critical superheating 
temperature and (b) melting temperature for different densities of 0.58, 0.78 and 1.13 g/cm^. 
Logarithm longitudinal scale is used. 

while only one at 1.13 g/cm^. In addition, for critical superheated solid, the peaks broaden 
while moving to higher energy as density is increased. All of the superheated solid sam- 
ples exhibit nearly free-electron characters. Again, the liquid spectrum is featureless and 
Drude-type. This allows us to expect the optical measurement to be able to diagnose the 
homogeneous melting. 

Through extrapolating to the zero frequency limit, the dc conductivity can be determined 
by fitting with the simple Drude fornt^ 



11 




12 - 







0.6 



0.7 



0.8 



0.9 



1.0 



1.1 



P (g/cm') 



Figure 8: (Color online) Electron dc conductivities at critical superheating temperatures (squares) 
and melting temperatures (circles) for different densities. 



where td represents effective collision time. Figure 8 presents the dc conductivity of lithium 
at Tis and for different densities. The dc conductivity shows a systematic behavior in 
terms of density for both cases. Stronger ion-electron scattering with increasing density 
would diminish the conductivity. In addition, the conductivity decreases with the temper- 
ature for metals. It is thus natural that the conductivity decrease with increasing density 
in the region with positive melting line slope. However, in the region with negative melting 
line slope, the effect of increased scattering is more prominent than that of the decreased 
temperature, which plays a central role in diminishing the conductivity. The dc conductivi- 
ties at Tm are in order of 10^ (ficm)"^, which is consistent with the experimental and other 
theoretical results of the liquid lithium in the similar density and temperature ranges^^*^. 
Also it can be seen from Fig. 8 that the dc conductivity difference between Tm and Tig 
becomes smaller with increasing density, and even undistinguished at 1.13g/cm^. 

IV. CONCLUSION 

In summary, we have performed FPMD simulations of the melting curve of lithium up to 
30 GPa using Z-method. The results are in good agreement with experimental measurements 
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and the two-phase simulations. It can be concluded that the melting line maximum of lithium 
may be caused by higher compressibihty of hquid phase than sohd phase, without LLPT 
accompanied. In addition, we have also systematically studied the atomic dynamic diffusion 
behavior and electronic dynamic conductivity properties at the critical superheating and 
melting points of lithium, which have revealed prominent physical differences between the 
superheated solid phase and the disordered liquid phase. For these two homogeneous phases, 
interestingly, the electron conductivities, especially the dc components, show the merging 
tendency at high densities, which suggests the increasing role the local structure plays in 
determining the electron-ion scattering. 
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